function [Y,Xx,Ea,Eu,ExitFlag] = run(S,IxPer)
% simulate.global.run  [Not a public function] Global nonlinear ...
% simulation in bkw models.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2015 IRIS Solutions Team.

%--------------------------------------------------------------------------

% Run first-order first
%-----------------------
S.M = simulate.linear.multipliers(S);
[S.y,S.x,S.Ea,S.Eu] = simulate.linear.run(S,1);

% Prepare global nonlinear
%--------------------------
ny = size(S.Z,1);
nx = size(S.T,1);
nb = size(S.T,2);
nf = nx - nb;
ne = size(S.Ea,1);
t = find(IxPer,1);
tt = 2;
ttAbs = -S.MinT + t;

% Unlog
%-------
isLog = ~isequal(S.Global.Unlog,@all);
if ~isLog
    S.x0(S.IxXLog(nf+1:end)) = exp( S.x0(S.IxXLog(nf+1:end)) );
    S.y(S.IxYLog) = exp( S.y(S.IxYLog) );
    S.x(S.IxXLog) = exp( S.x(S.IxXLog) );
end

eqtnN = S.Global.EqtnN;
ixXCurr = S.Global.IxXCurr;
g = S.Global.Times;

pp = real(S.Assign(1,S.NameType==4));
LL = S.L;

ixY = ~S.Anch(1:ny,1);
ixX = ~S.Anch(ny+(1:nx),1) & ixXCurr;
ixEa = S.Anch(ny+nx+(1:ne),1);
ixEu = S.Anch(ny+nx+ne+(1:ne),1);

nzy = sum(ixY);
nzx = sum(ixX);
nzea = sum(ixEa);

% Init condition from first-order simulation.
z0 = [S.y(ixY,1);S.x(ixX,1);S.Ea(ixEa,1);S.Eu(ixEu,1)];

% Call Optim Tbx
%----------------
if isequal(S.Solver,@fsolve)
    [z1,res,ExitFlag] = ...
        fsolve(@doObjFunc,z0,S.OptimSet);
elseif isequal(S.Solver,@lsqnonlin)
    [z1,~,res,ExitFlag] = ...
        lsqnonlin(@doObjFunc,z0,[],[],S.OptimSet);
end
if ExitFlag < 0 || any(abs(res) > S.Global.Tolerance)
    ExitFlag = -10;
end

if ExitFlag <= 0
    doFail();
end

% Output data
%-------------
if any(ixY)
    S.y(ixY,1) = z1(1:nzy);
end
if any(ixX)
    S.x(ixX,1) = z1(nzy+(1:nzx));
end
if any(ixEa)
    S.Ea(ixEa,1) = z1(nzy+nzx+(1:nzea));
end
if any(ixEu)
    S.Eu(ixEu,1) = z1(nzy+nzx+nzea+1:end);
end

if ~isLog
    S.x0(S.IxXLog(nf+1:end)) = log( S.x0(S.IxXLog(nf+1:end)) );
    S.y(S.IxYLog) = log( S.y(S.IxYLog) );
    S.x(S.IxXLog) = log( S.x(S.IxXLog) );
end

Y = S.y;
Xx = S.x;
Ea = S.Ea;
Eu = S.Eu;

return

    
    
    function D = doObjFunc(Z)
        % Measurement variables.
        y1 = [ nan(ny,1), S.y(:,1) ];
        if any(ixY)
            y1(ixY,2) = Z(1:nzy);
        end
        
        % Transition variables.
        xx1 = [ [nan(nf,1);S.x0], S.x(:,1) ];
        if any(ixX)
            xx1(ixX,2) = Z(nzy+(1:nzx));
        end
        
        if S.IsDeviation && S.IsAddSstate
            if ~isempty(y1)
                y1 = y1 + S.YBar(:,1+(t-1:t));
            end
            xx1 = xx1 + S.XBar(:,1+(t-1:t));
        end
        
        if isLog
            if any(S.IxYLog)
                y1(S.IxYLog,:) = real(exp(y1(S.IxYLog,:)));
            end
            if any(S.IxXLog)
                xx1(S.IxXLog,:) = real(exp(xx1(S.IxXLog,:)));
            end
        end
        
        % Anticipated shocks.
        ea1 = [ nan(ne,1), S.Ea(:,1) ];
        if any(ixEa)
            ea1(ixEa,2) = Z(nzy+nzx+(1:nzea));
        end
        
        % Unanticipated shocks.
        eu1 = [ nan(ne,1), S.Eu(:,1) ];
        if any(ixEu)
            eu1(ixEu,2) = Z(nzy+nzx+nzea+1:end);
        end
        e1 = ea1 + eu1;
        
        try
            % Evaluate all equations at once
            %--------------------------------
            D = g*eqtnN(y1,xx1,e1,pp,tt,LL,ttAbs);
        catch
            % Evaluate individual equations to catch errors
            %----------------------------------------------
            nEqtn = length(S.Eqtn);
            errMsg = { };
            for i = 1 : nEqtn
                if isempty(S.Global.EqtnNI{i})
                    continue
                end
                try
                    S.Global.EqtnNI{i}(y1,xx1,e1,pp,tt,LL,ttAbs);
                catch Err
                    errMsg = [errMsg, {t,S.Eqtn{i},Err.message} ]; %#ok<AGROW>
                end
            end
            if ~isempty(errMsg)
                utils.error('simulate:global:run', ...
                    ['Error evaluating equations when ', ...
                    'running global nonlinear simulation at t=%g ', ...
                    'in this equation: %s\n', ...
                    '\tUncle says: %s'], ...
                    errMsg{:});
            end
        end
    end % doObjFuncXQ()



    function doFail()
        if S.IsError
            % @@@@@ MOSW
            msgFunc = @(varargin) utils.error(varargin{:});
        else
            % @@@@@ MOSW
            msgFunc = @(varargin) utils.warning(varargin{:});
        end
        [~,pos] = max(abs(res));
        pos = pos(1);
        msgFunc('simulate:global:run', ...
            ['Global nonlinear simulation ', ...
            'may be inaccurate at t=%g. ', ...
            'Max abs discrepancy %g in this equation: %s'], ...
            t,max(abs(res)),S.Eqtn{pos});
    end % doFail()
end % main
